
******************************************
*     Code for Ahlskog (2023), JEPS      *
* Replication of Rasmussen et al. (2021) *
******************************************

/*
Important variables:
LopNrParID is the twin pair identifier
TwinNr identifies twin nr 1/2 within a pair (arbitrary)
zyg is the zygosity, with 1 being MZ and 2 being DZ
eduYears is education years
PGI_EA_single is a polygenic index of education

*/

set scheme s1mono
cd "C:\Userdata\Shared\Dofiles\PrelDoFiles\PrelRafael\Rasmussen replication\Output"
log using "logfile.txt", replace
use "C:\Userdata\Shared\Dofiles\PrelDoFiles\PrelRafael\Rasmussen replication\main_dataset", clear








****************************
* CREATE IDEOLOGY OUTCOMES *
****************************


* Economic

* Scale based on original
gen original_econlib = decrease_public_sector + decrease_welfare + more_freedom_companies + lower_taxes - decrease_economic_inequality


* Rescale variables negative related to econlib
* (makes the calculation of the constraint measure easier below)
replace decrease_finmarket_impact = 6 - decrease_finmarket_impact
replace keep_property_taxes = 6 - keep_property_taxes
replace keep_maxtaxa = 6 - keep_maxtaxa
replace decrease_economic_inequality = 6 - decrease_economic_inequality

* First, conduct sem of all the suspected econlib variables in the MZ data
local ELVARS = "decrease_public_sector decrease_welfare lower_taxes sell_public_enterprise more_private_healthcare more_free_trade more_freeschools more_freedom_companies decrease_finmarket_impact keep_property_taxes keep_maxtaxa decrease_economic_inequality"
sem (`ELVARS'  <- EL) if zyg == 2 & TwinNr == 1, method(ml)
matrix ELWEIGHTS = e(b)

* Then, use the coefficients from the sem to construct the econlib outcome, and rescale back to 1-5
gen econlib = 0
local i = 1
foreach v in `ELVARS' {
	replace econlib = econlib + ELWEIGHTS[1,`i'] * `v'
	local i = `i' + 2
}
su econlib
replace econlib = 4 * ((econlib - r(min)) / (r(max) - r(min))) + 1



* Do the same but in DZ twins
sem (`ELVARS'  <- EL) if zyg == 1 & TwinNr == 1, method(ml)
matrix ELWEIGHTS = e(b)

* Then, use the coefficients from the sem to construct the econlib outcome, and rescale back to 1-5
gen econlib_DZ = 0
local i = 1
foreach v in `ELVARS' {
	replace econlib_DZ = econlib_DZ + ELWEIGHTS[1,`i'] * `v'
	local i = `i' + 2
}
su econlib_DZ
replace econlib_DZ = 4 * ((econlib_DZ - r(min)) / (r(max) - r(min))) + 1




* Social

* Scale based on original
gen original_conserv = harder_punishment + language_test_citizenship - more_support_immigrant_culture - less_carbondioxide - decrease_pollution


* Rescale for convenience
replace more_skilled_immigration = 6 - more_skilled_immigration
replace more_support_immigrant_culture = 6 - more_support_immigrant_culture
replace decrease_defense_spending = 6 - decrease_defense_spending

* First, conduct sem of all the suspected econlib variables in the MZ data
local CONVARS = "earlier_grades ban_pornography limit_abortion harder_punishment language_test_citizenship fewer_refugees more_skilled_immigration more_support_immigrant_culture decrease_defense_spending"
sem (`CONVARS'  <- CON) if zyg == 2 & TwinNr == 1, method(ml)
matrix CONWEIGHTS = e(b)

* Then, use the coefficients from the sem to construct the econlib outcome, and rescale back to 1-5
gen conserv = 0
local i = 1
foreach v in `CONVARS' {
	replace conserv = conserv + CONWEIGHTS[1,`i'] * `v'
	local i = `i' + 2
}
su conserv
replace conserv = 4 * ((conserv - r(min)) / (r(max) - r(min))) + 1


* Same but in DZ twins
sem (`CONVARS'  <- CON) if zyg == 1 & TwinNr == 1, method(ml)
matrix CONWEIGHTS = e(b)

* Then, use the coefficients from the sem to construct the econlib outcome, and rescale back to 1-5
gen conserv_DZ = 0
local i = 1
foreach v in `CONVARS' {
	replace conserv_DZ = conserv_DZ + CONWEIGHTS[1,`i'] * `v'
	local i = `i' + 2
}
su conserv_DZ
replace conserv_DZ = 4 * ((conserv_DZ - r(min)) / (r(max) - r(min))) + 1






**************
* OTHER PREP *
**************

* Labels
label variable PGI_EA_single "PGI EA"
label variable eduYears "Education years"
label variable original_econlib "Orig. EI"
label variable original_conserv "Orig. SI"
label variable econlib "New EI"
label variable conserv "New SI"
label variable econlib_DZ "New EI"
label variable conserv_DZ "New SI"

* Standardize PGI by birth year
foreach v in PGI_EA_single {
	egen temp_mean = mean(`v'), by(BirthYear)
	egen temp_sd = sd(`v'), by(BirthYear)
	replace `v' = (`v'-temp_mean)/temp_sd
	drop temp_mean temp_sd
}


preserve
	keep if zyg==1
	keep if original_conserv!=. | original_econlib!=.
	keep BirthYear sex eduYears original_econlib original_conserv econlib conserv earlier_grades ban_pornography limit_abortion harder_punishment language_test_citizenship fewer_refugees more_skilled_immigration more_support_immigrant_culture decrease_defense_spending decrease_public_sector decrease_welfare lower_taxes sell_public_enterprise more_private_healthcare more_free_trade more_freeschools more_freedom_companies decrease_finmarket_impact keep_property_taxes keep_maxtaxa decrease_economic_inequality less_carbondioxide decrease_pollution
	order BirthYear sex eduYears original_econlib original_conserv econlib conserv earlier_grades ban_pornography limit_abortion harder_punishment language_test_citizenship fewer_refugees more_skilled_immigration more_support_immigrant_culture decrease_defense_spending decrease_public_sector decrease_welfare lower_taxes sell_public_enterprise more_private_healthcare more_free_trade more_freeschools more_freedom_companies decrease_finmarket_impact keep_property_taxes keep_maxtaxa decrease_economic_inequality less_carbondioxide decrease_pollution
	outreg2 using "Table_A2.text", replace title("Descriptive statistics, MZ") tex(frag) label sum(log)
restore

preserve
	keep if zyg==2
	keep if ( original_conserv!=. | original_econlib!=. ) & PGI_EA_single!=.
	keep BirthYear sex eduYears original_econlib original_conserv econlib_DZ conserv_DZ PGI_EA_single earlier_grades ban_pornography limit_abortion harder_punishment language_test_citizenship fewer_refugees more_skilled_immigration more_support_immigrant_culture decrease_defense_spending decrease_public_sector decrease_welfare lower_taxes sell_public_enterprise more_private_healthcare more_free_trade more_freeschools more_freedom_companies decrease_finmarket_impact keep_property_taxes keep_maxtaxa decrease_economic_inequality less_carbondioxide decrease_pollution
	order BirthYear sex eduYears original_econlib original_conserv econlib_DZ conserv_DZ PGI_EA_single earlier_grades ban_pornography limit_abortion harder_punishment language_test_citizenship fewer_refugees more_skilled_immigration more_support_immigrant_culture decrease_defense_spending decrease_public_sector decrease_welfare lower_taxes sell_public_enterprise more_private_healthcare more_free_trade more_freeschools more_freedom_companies decrease_finmarket_impact keep_property_taxes keep_maxtaxa decrease_economic_inequality less_carbondioxide decrease_pollution
	outreg2 using "Table_A3.text", replace title("Descriptive statistics, DZ") tex(frag) label sum(log)
restore



* Standardize variables and keep only complete pairs
foreach v in earlier_grades ban_pornography limit_abortion harder_punishment ///
	language_test_citizenship fewer_refugees more_skilled_immigration ///
	more_support_immigrant_culture decrease_defense_spending decrease_public_sector ///
	decrease_welfare more_freedom_companies lower_taxes decrease_economic_inequality ///
	eduYears econlib conserv original_econlib original_conserv econlib_DZ conserv_DZ {
	
	su `v'
	replace `v' = (`v' - r(mean)) / r(sd)
	bysort LopNrParID: replace `v' = . if `v'[1] == . | `v'[2] == .
}








************
* ANALYSES *
************


* MZ models - corresponds to table 1 in the main paper
reghdfe original_econlib eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_1.tex", replace label tex(frag) addtext(Twin pair FE, YES)
reghdfe econlib eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_1.tex", append label tex(frag) addtext(Twin pair FE, YES)
reghdfe original_conserv eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_1.tex", append label tex(frag) addtext(Twin pair FE, YES)
reghdfe conserv eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_1.tex", append label tex(frag) addtext(Twin pair FE, YES)


* DZ models - corresponds to table 3 in the main paper
reghdfe original_econlib PGI_EA_single sex if zyg == 2, cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_3.tex", replace label tex(frag) addtext(Twin pair FE, YES)
reghdfe econlib_DZ PGI_EA_single sex if zyg == 2, cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_3.tex", append label tex(frag) addtext(Twin pair FE, YES)
reghdfe original_conserv PGI_EA_single sex if zyg == 2, cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_3.tex", append label tex(frag) addtext(Twin pair FE, YES)
reghdfe conserv_DZ PGI_EA_single sex if zyg == 2, cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_3.tex", append label tex(frag) addtext(Twin pair FE, YES)



* What drives differences? - corresponds to table 2 in the main paper
file open itemtable using "Table_2.tex", write replace
file write itemtable "\begin{tabular}{lcc} \\" _n
file write itemtable "\multicolumn{3}{c}{\textbf{Approximation of original items on Economic Ideology}} \\ \hline" _n

	* This shows that results for original_econlib are entirely driven by lower_taxes
	reghdfe decrease_public_sector eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Decrease the public sector & `b' & (`se') \\" _n
	reghdfe decrease_welfare eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Decrease social welfare & `b' & (`se') \\" _n
	reghdfe more_freedom_companies eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Give companies more freedom & `b' & (`se') \\" _n
	reghdfe lower_taxes eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Taxes should be cut & `b'** & (`se') \\" _n
	reghdfe decrease_economic_inequality eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Decrease economic inequality (reversed) & `b' & (`se') \\ \\" _n _n

	* ...while results for the new conserv scale are driven by migration related attitudes
	reghdfe earlier_grades eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Give school grades at younger age & `b' & (`se') \\" _n
	reghdfe ban_pornography eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Ban pornography & `b' & (`se') \\" _n
	reghdfe limit_abortion eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Limit access to abortion & `b' & (`se') \\" _n
	reghdfe harder_punishment eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Introduce much harder punishment for criminals & `b'* & (`se') \\" _n
	reghdfe language_test_citizenship eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Instate a language test for Swedish citizenship & `b'* & (`se') \\" _n
	reghdfe fewer_refugees eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Admit fewer refugees & `b'*** & (`se') \\" _n
	reghdfe more_skilled_immigration eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Allow more skilled immigration (reversed) & `b'*** & (`se') \\" _n
	reghdfe more_support_immigrant_culture eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Increase assistance for immigrants' culture (reversed) & `b'** & (`se') \\" _n
	reghdfe decrease_defense_spending eduYears if zyg == 1, cluster(LopNrParID) absorb(LopNrParID)
		local b : di %12.4f _b[eduYears]
		local b = strtrim("`b'")
		local se : di %12.4f _se[eduYears]
		local se = strtrim("`se'")
		file write itemtable "Decrease defense spending (reversed) & `b' & (`se') \\ \\" _n _n
	
	file write itemtable "\multicolumn{3}{p{0.8\textwidth}}{Notes: All variables are standardized. Robust standard errors in parentheses. *** p$<$0.01, ** p$<$0.05, * p$<$0.1} \\"
	file write itemtable "\end{tabular}"
	file close itemtable




* Main models with matching sample sizes
* MZ models - corresponds to table A4 in the online appendix
reghdfe original_econlib eduYears if zyg == 1 & original_econlib!=. & econlib!=. & original_conserv!=. & conserv!=. , cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_A4.tex", replace label tex(frag) addtext(Twin pair FE, YES)
reghdfe econlib eduYears if zyg == 1 & original_econlib!=. & econlib!=. & original_conserv!=. & conserv!=., cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_A4.tex", append label tex(frag) addtext(Twin pair FE, YES)
reghdfe original_conserv eduYears if zyg == 1 & original_econlib!=. & econlib!=. & original_conserv!=. & conserv!=., cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_A4.tex", append label tex(frag) addtext(Twin pair FE, YES)
reghdfe conserv eduYears if zyg == 1 & original_econlib!=. & econlib!=. & original_conserv!=. & conserv!=., cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_A4.tex", append label tex(frag) addtext(Twin pair FE, YES)


* DZ models - corresponds to table A5 in the online appendix
reghdfe original_econlib PGI_EA_single sex if zyg == 2 & original_econlib!=. & econlib!=. & original_conserv!=. & conserv!=., cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_A5.tex", replace label tex(frag) addtext(Twin pair FE, YES)
reghdfe econlib_DZ PGI_EA_single sex if zyg == 2 & original_econlib!=. & econlib!=. & original_conserv!=. & conserv!=., cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_A5.tex", append label tex(frag) addtext(Twin pair FE, YES)
reghdfe original_conserv PGI_EA_single sex if zyg == 2 & original_econlib!=. & econlib!=. & original_conserv!=. & conserv!=., cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_A5.tex", append label tex(frag) addtext(Twin pair FE, YES)
reghdfe conserv_DZ PGI_EA_single sex if zyg == 2 & original_econlib!=. & econlib!=. & original_conserv!=. & conserv!=., cluster(LopNrParID) absorb(LopNrParID)
outreg2 using "Table_A5.tex", append label tex(frag) addtext(Twin pair FE, YES)




* Prep R dataset for calculating Bayes factors

preserve
	bysort LopNrParID: gen diff_edu = eduYears[1] - eduYears[2]
	replace diff_edu=. if BESTZYG!=1
	bysort LopNrParID: gen diff_PGI = PGI_EA_single[1]-PGI_EA_single[2]
	replace diff_PGI=. if BESTZYG==1
	bysort LopNrParID: gen diff_econlib = econlib[1]-econlib[2]
	bysort LopNrParID: gen diff_conserv = conserv[1]-conserv[2]
	bysort LopNrParID: gen diff_econlib_DZ = econlib_DZ[1]-econlib_DZ[2]
	bysort LopNrParID: gen diff_conserv_DZ = conserv_DZ[1]-conserv_DZ[2]
	bysort LopNrParID: gen diff_orig_econlib = original_econlib[1]-original_econlib[2]
	bysort LopNrParID: gen diff_orig_conserv = original_conserv[1]-original_conserv[2]
	keep if TwinNr==1
	keep diff*
	save "C:\Userdata\Shared\Dofiles\PrelDoFiles\PrelRafael\Rasmussen replication\bfdata.dta", replace
restore



log close
